#######################################
######Masculinity Threat Analysis######
#######################################
###############10/20/25##################
#######################################

####Analysis Script for Gothreau & Haas, forthcoming in Journal of Experimental Political Science###
####A Replication and Extension of Willer et al. (2013) Overdoing Gender: A Test of the Masculine Overcompensation Thesis###

library(haven)
library(ggplot2)
library(dplyr)
library(broom)
library(officer)
library(flextable)
library(forcats)
library(psych)
library(survey)
library(stargazer)
library(purrr)
library(tidyr)
library(forcats)

#Set Working Directory#
MTData <- read.csv("CleanMTData.csv")

# Recode factors -------------------------------------------------------------
MTData$P_EXP <- factor(
  MTData$P_EXP,
  levels = c("threat","no threat","constrained","PCKT"),
  labels = c("threat","no threat","constrained","PCKT")
)

MTData$gender_binary <- factor(
  MTData$Q_CDGDRDESC_17,
  levels = c("1","2"),
  labels = c("Men","Women")
)

# Create datasets by gender
Data_female <- subset(MTData, gender_binary == "Women")
Data_male   <- subset(MTData, gender_binary == "Men")


##T-test for additional cancel culture item (all conditions except control received this item)

dvs <- "CANCEL_CULTURE_rescaled"

run_weighted_ttests <- function(data, label, weight_var = "WEIGHT") {
  # compare each treatment to PCKT; EXCLUDE "no threat"
  experimental_groups <- setdiff(unique(data$P_EXP), c("PCKT", "no threat"))
  results <- list()
  
  for (dv in dvs) {
    # ensure DV is numeric
    if (!dv %in% names(data)) next
    if (!is.numeric(data[[dv]])) {
      suppressWarnings({ data[[dv]] <- as.numeric(data[[dv]]) })
    }
    
    for (group in experimental_groups) {
      subset_data <- data %>%
        filter(P_EXP %in% c("PCKT", group)) %>%
        select(P_EXP, !!sym(weight_var), !!sym(dv)) %>%
        drop_na()
      
      # must have rows from both groups
      subset_data$P_EXP <- droplevels(subset_data$P_EXP)
      if (nlevels(subset_data$P_EXP) != 2) next
      
      # set reference and build design
      subset_data$P_EXP <- relevel(subset_data$P_EXP, ref = "PCKT")
      design <- svydesign(ids = ~1, weights = as.formula(paste0("~", weight_var)), data = subset_data)
      
      form <- as.formula(paste(dv, "~ P_EXP"))
      ttest <- try(svyttest(form, design), silent = TRUE)
      if (inherits(ttest, "try-error")) next
      
      # weighted means
      group_means <- svyby(as.formula(paste0("~", dv)), ~P_EXP, design, svymean, na.rm = TRUE)
      mean_PCKT  <- as.numeric(group_means[group_means$P_EXP == "PCKT", 2])
      mean_treat <- as.numeric(group_means[group_means$P_EXP == group, 2])
      
      out <- data.frame(
        gender = label,
        DV = dv,
        experimental_group = group,
        mean_PCKT = round(mean_PCKT, 2),
        mean_treat = round(mean_treat, 2),
        t = round(unname(ttest$statistic), 2),
        p = unname(ttest$p.value),
        sig = ifelse(unname(ttest$p.value) < .001, "***",
                     ifelse(unname(ttest$p.value) < .01,  "**",
                            ifelse(unname(ttest$p.value) < .05,  "*", ""))),
        stringsAsFactors = FALSE
      )
      
      results[[paste(label, dv, group, sep = "_vs_")]] <- out
    }
  }
  
  dplyr::bind_rows(results)
}

ttests_female <- run_weighted_ttests(Data_female, "female")
ttests_male   <- run_weighted_ttests(Data_male,   "male")
all_ttests    <- dplyr::bind_rows(ttests_female, ttests_male)
all_ttests

###Table for Appendix

# Order groups as you prefer in the table
group_order <- c("threat", "constrained", "PCKT")

table_all <- all_ttests %>%
  filter(experimental_group != "no threat") %>%
  mutate(
    DV = factor(DV, levels = dvs),
    experimental_group = factor(experimental_group, levels = group_order),
    mean_treat = round(mean_treat, 2),
    comparison = paste0("t = ", round(t, 2), sig, ", p = ", sprintf("%.3f", p))
  ) %>%
  select(gender, DV, experimental_group, mean_treat, comparison) %>%
  arrange(gender, DV, experimental_group)

# ---- Pretty flextable with all groups ----
ft <- flextable(table_all)
ft <- set_header_labels(
  ft,
  gender            = "gender",
  DV                = "DV",
  experimental_group= "Group",
  mean_treat        = "Mean (Group)",
  comparison        = "t-test vs. PCKT"
)
ft <- theme_vanilla(ft)
ft <- flextable::align(ft, j = c("gender", "DV", "experimental_group"), align = "left")
ft <- flextable::align(ft, j = c("mean_treat", "comparison"), align = "center")
ft <- autofit(ft)

# Optional: visually separate gender blocks
# (adds a light line where gender changes)
if (nrow(table_all) > 1) {
  brks <- which(diff(as.integer(factor(table_all$gender))) != 0)
  if (length(brks)) {
    ft <- hline(ft, i = brks, border = fp_border(color = "#CCCCCC"), part = "body")
  }
}

# Legend for significance stars
ft <- add_footer_lines(ft, "* p < .05, ** p < .01, *** p < .001")

# ---- Export to Word ----
doc <- read_docx()
doc <- body_add_par(doc, "T-Test Results: All Experimental Groups vs. PCKT", style = "heading 1")
doc <- body_add_flextable(doc, ft)
print(doc, target = "Study1_replication_all_groups_vs_PCKT.docx")
